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Engineer  Topographic  Laboratories,  Fort  Belvoir,  Virginia  under 
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Project  Engineer  for  the  government. 
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Introduction 


This  report  describes  the  mathematical  analysis  on  which 
the  program  modifications  to  Fotonap  are  based.  Implementation  of 
the  option  not  to  compute  the  covariance  matrix  of  the  solution 
vector  did  not  involve  any  new  analysis  and  therefore  is  not 
described  here.  It  did,  however,  involve  the  programming  of  a 
completely  new  back-substitution  subroutine  (INVRTC) . The  time 
saving  achieved  when  this  option  was  exercised  was  also  quite 
significant,  the  time  taken  to  obtain  the  solution  of  the  normal 
equations  matrix  being  reduced  from  1 hour  50  minutes  to  20  minutes 
for  a typical  run  using  photogrammetric  measurements.  The  main 
part  (Section  2)  of  this  report  describes  Geoceiver  measurements. 

Section  3 deals  with  the  Hopfield  tropospheric  refraction 
correction  formula.  Although  the  formula  used  in  Fotonap  is 
algebraically  similar  to  the  NSWC  version  of  the  Hopfield  model,  it 
is  computationally  rather  different.  Also,  in  the  Fotonap  version 
of  the  Hopfield  model,  default  values  for  temperature  and  pressure 
are  adjusted  for  station  height. 

The  Fotonap  User's  Guide  has  been  modified  to  be  consistent 
with  the  new  version  of  the  program.  The  changes  to  the  User's 
Guide  are  described  in  Section  4. 

As  the  changes  to  Fotonap  to  handle  Geoceiver  measurements 
had  to  be  fairly  extensive,  some  changes  to  the  program  structure  for 
computing  predicted  measurements  (involving  subroutines  PRTIAL, 
MEASUR,  MESOLD)  were  made  in  order  to  facilitate  such  changes.  In 
order  to  check  that  these  changes  did  not  introduce  errors  in  the 
computation  of  existing  measurement  types,  a test  case  involving  all 


in 


different  measurement  types  was  devised.  Interestingly,  this  revealed 
errors  in  the  computation  of  some  of  the  partial  derivatives  of 
the  old  version  of  the  program.  The  analysis  for  these  changes 
to  program  are  not  included  in  this  report,  but  are  based  on  the 
analysis  in  reference  9. 


2. 


Geoceiver  Measurements 


Although  geoceiver  measurements  may  he  regarded  as 
measurements  of  range  differences,  the  correspondence  between  the 
two  is  exact  only  in  the  absence  of  an  atmosphere  between  the 
transmitter  and  the  receiver.  In  practice,  of  course,  the  atmos- 
phere cannot  be  removed.  However,  the  atmospheric  effect  may  be 
removed.  The  two  parts  of  the  atmosphere,  the  ionosphere  and  the 
troposphere,  have  quite  different  effects  on  the  measurements  and 
must  therefore  be  treated  separately.  The  ionospheric  effect, 
which  is  frequency  dependent,  may  be  estimated  through  the  use  of 
a two-frequency  transmission.  The  tropospheric  effect  may  be 
cn  ated  with  the  aid  of  a tropospheric  model.  In  order  to  do 
wever,  it  is  necessary  to  know  the  geometrical  relation- 
owips  between  the  transmitter,  the  receiver  and  the  troposphere. 
For  that  reason,  tropospheric  corrections  are  usually  computed 
in  an  orbit  determination  program  rather  than  in  a preprocessor. 
(Section  3 of  this  report  describes  the  NWL-Hopfield  tropospheric 
model  and  how  it  is  used  in  Photonap  to  compute  tropospheric 
corrections.)  In  the  following  description  of  geoceiver  measure- 
ments it  will  be  asstimed  that  atmospheric  corrections  have  been 
made  to  the  data.  For  this  reason  no  further  mention  will  be 
made  of  the  atmosphere. 

A satellite  transmits  a constant  frequency  signal  that 
is  received  by  the  geoceiver.  The  received  frequency  will, 
because  of  the  motion  of  the  transmitter  relative  to  the  receiver, 
differ  from  the  transmission  frequency  by  the  (one-way)  Doppler 
frequency. 
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Let  R(t)  denote  the  retarded  range  of  the  transmitter 
relative  to  receiver  at  time  t.  In  other  words,  R(t)  is  the 
length  of  the  signal  path  for  a signal  being  received  at  time  t. 
Furthermore,  let  R(t)  denote  the  derivative  of  R(t)  with  respect 
to  t,  and  let  c denote  the  transmission  velocity.  Then  clearly 


(2.1) 


where  Vo  is  the  transmission  frequency,  and  6v(t)  is  the  Doppler 
frequency  at  time  t.  Integrating  the  above  equation  between  times 
(t  - T)  and  t,  we  deduce  that 


g(t)  = R(t)  - R(t-T),  (2.2) 

where  the  geoceiver  measurement  g(t)  is  given  by 

g(t)  = ^ 6v(t  + t - T)  dx.  (2.3) 

•'o 

The  geoceiver  equipment  may  be  either  ground  based  or  carried  by 
a satellite.  The  first  case  is  the  more  common  one  and  will  be 
considered  first. 

2. 1 Satellite- to-Ground  Measurements 

Let  Xso (t)  denote  the  satellite  position  vector  at  time 
t in  the  'Mean  of  date,  1950.0'  coordinate  system.  (This  is  the 
inertial  coordinate  system  used  by  Photonap  for  all  orbit  computa- 
tions.) Also,  let  Rd5  denote  the  rotation  matrix  transforming 
vectors  in  the  '1950.0'  system  to  the  'true  of  date  system'  at 
time  t.  Although  Ros  obviously  is  a function  of  time,  it  may  be 
regarded  as  constant  over  a period  of  about  a minute.  The  rotation 
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matrix  Rfo(t) , which  transforms  an  'of  date'  vector  to  an  'Earth 
fixed'  vector  is  defined  by 


Rfo(t)  = 


cosA(t) 

-sinA(t) 

0 


sinA(t)  0" 
cosA(t)  0 

0 ij. 


(2.4) 


where  A(t)  is  the  Greenwich  nour  angle,  and  A(t)  = w , the  Earth's 
rotation  rate. 

The  inertial  satellite  vector  x(x,t)  is  defined  by 


x(T,t)  = Rfo(t)  Ro5Xso(t) 


(2.5) 


Equation  (2.5)  can  be  seen  to  represent  the  satellite  vector  in 
an  inertial  coordinate  system  that  coincides  with  the  'Earth 
fixed'  coordinate  system  at  time  t = t. 

Denoting  the  'Earth  fixed'  station-vector  by  Sf , the 
retarded  range,  R(t) , must,  by  definition,  satisfy  the  equation 

R(t)  =\/[s,  - 5(T,t)]  • [sf  - x(T,t)]  , (2.6) 

where,  t = t - R(t)/c  (2.7) 

The  simultaneous  equations  (2.6)  and  (2.7)  cannot  be  solved 
explicitly.  However,  on  account  of  R(t)/c  being  small,  the 
equations  can  very  easily  be  solved  using  an  iterative  technique. 

The  predicted  geoceiver  measurement  is  then  computed 
using  equation  (2.2). 
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2.1.1  Measurement  Partial  Derivatives 
Let  q(t)  be  defined  by 

q(t)  = X5o(t)  (2.8) 

where  p is  any  parameter  affecting  the  orbit.  (p  may,  e.g. , 
be  a component  of  the  initial  position  vector  or  a coefficient 
of  one  of  the  spherical  harmonics  of  the  gras'ity  field.  Photonap 
computes  a vector  q for  each  unknown  parameter  affecting  the  orbit.) 

Squaring  both  sides  of  equation  (2.6)  and  differentiating 
with  respect  to  p we  obtain 

R(t)  R(t)  - S,]-[|^  J(T,t)  + If  |1]  (2.9) 

Writing 

u =[x(T.t)  - Sp]/R(t)  (2.10) 

and 

Ro  = u • Iy  x(T,t)  (2.11) 

we  find  with  the  aid  of  equation  (2.7)  that 

Ip  ««  = S • Ip  *(^.«  + R"!-  I 1^  K(t)]  (2.12) 

It  hence  follows  from  the  above  and  equations  (2.5)  and  (2.8)  that 
R(t)  = Uo  ■ Rfo(t)  Rosq  (t),  (2.13) 


where 

Uo  = u/ (1  + Ro/c) . 


(2.14) 


I 
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Differentiating  the  square  of  equation  (2.6)  with  respect 
to  the  station  vector  we  similarly  obtain 


R(t)  -|g  R(t)  -j^x(T,t)  - x(T,t)  ’ 


whence 


ls,R('>  ■ R"  Mlf,R<‘>]  - "• 


We  hence  deduce  that 


3S. 


R(t)  = - Uc 


(2.15) 


To  obtain  the  time  derivative  R(t) , we  similarly  deduce 
from  equation  (2.6)  that 


R(t)R(t)  = Jx(x,t)  - Sj] ■ [ It  It  x(T,t)  llj ^ 


whence  by  equations  (2.10),  (2.11)  and  (2.7) 

R(t)  = ^ • It  x(T,t)  + Ro^l  - ^ R(t)^ 
From  the  above  equation  it  follows  that 


R(t)  = Ro/(l  + Ro/c)  + Uo  • It  x(x,t) 


(2.16) 


It  may  readily  be  verified  that  with  the  rotation  matrix  R^^^  given 
by  equation  (2.4) 


R..  R?-  = u) 


fO  FD 


0 10 

-10  0 
0 0 0 


(2.17) 


Since  it  follows  from  equation  (2.5)  that 


^ i(x,t)  = R,^(t)  R^^  X5o(x). 
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and 


Ro5  X50  (t)  = (t)  X(T,t), 


we  deduce  with  the  aid  of  equation  (2.17)  that 


It 


X2 (T ,t) 
■Xi (x ,t) 
0 


Furthermore,  since 


(Xi-Si)X2  - (X2-S2>Xl  = (Xi-Sl)S2  - (X2-S2)Si, 

it  follows  from  equation  (2.16)  that 

R(t)  = Ro/(l  + Ro/c)  + w|^(uo)i  S2  - (uo)2  SiJ 


(2.18) 


(2.19) 


2.2 


Satellite- to-Satellite  Measurements 


In  this  type  of  measurement  the  ground  station  is  replaced 
by  the  satellite  and  the  satellite  is  replaced  by  another  satelii’-e, 
whose  orbit  in  general  can  be  considered  known.  In  this  case, 
equation  (2.6)  is  replaced  by 


R(t)=-^[^X50  (t)  - Uso  (t)]  • [x50  (t)  " U50  (x)] 

where  Uso (x ) denotes  the  position  vector  of  the  transmitting 
satellite,  and  x satisfies  the  equation 

T = t - R(t)/c 


(2.20) 


(2.21) 


2.2.1  Measurement  Partial  Derivatives 

Similarly  to  the  satellite- to-ground  case,  we  find  that 
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R(t)  ^ R(t) 


Hence,  writing 


= 1^X50  (t)  - uso  (T)j  • jq(t)+  Gso  I 


V = [uso  (t)  - X50  (t)J/R(t), 


and 


Rj  = V • Uso 


we  deduce  that 


where 


R(t)  = -Vo  • q (t) 


Vo  = v/ (I  + R./c) 


Corresponding  to  equation  (2.8),  we  define 
q (t)  = |-  U50  (t)  . 

S >15 

It  then  follows  that 

. -V  - G„(t)  (- 

i.e.  , 

^ . q^(x)  - R^(l|p^R(t)). 

Hence , 

R(t)  - V.  • (t) 

I^Note  that  if  both  satellites  are  functions  of  the  same 
such  as  a gravity  coefficient,  so  that  p and  p^  are  the 

meter,  then  the  partial  derivative  is  simply  the  sum  of 
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(2.22) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 


(2.27) 

parameter, 
same  para- 
the  right 


hand  sides  of  equations  (2.24)  and  (2.27)j 

Differentiating  equation  (2.20),  we  obtain  for  the  time 


derivative 


R(t)  = - V 


1^x50  - u 


50  (t)(1  - J R(t)) 


i.e.  , 


R(t)  = - V • X50  + 


R,[l  - 


Consequently, 


R(t)  = - Vo  • X50  + Re /(I  + 


(2.28) 


3.0  The  NWL-Hopfleld  Model  for  Tropospheric  Refraction  of 
Radio  Ranging  Signals 

3. 1 The  Index  of  Refraction 

The  NWL-Hopfield  model,  like  many  similar  models,  is 
based  on  the  Smith-Weintraub  formula  (reference  1) 

N = [p  + 4810  l],  (3.1) 

where 

N is  the  refractivity  {N  = (n-l)lO®,  where  n is  the 
refractive  index} , 

T is  the  temperature  (degrees  Kelvin) , 
p is  the  total  pressure  (millibars) , 

e is  the  partial  pressure  of  water  vapor  (millibars) . 
Smith  and  Weintraub  claim  that  for  frequencies  below  30,000  MHz 
the  formula  is  accurate  to  within  0.5%  for  temperatures  between 
-50°C  and  40°C,  pressures  between  200  mb  and  1100  mb,  and  partial 
water  vapor  pressures  between  0 mb  and  30  mb. 

It  is  interesting  to  note  that  the  formula  is  derived 
from  the  three-term  expression 

N = ki  ^ + kz  I + ks  I2,  (3.2) 

where,  "the  first  term  expresses  the  sums  of  the  distortions  of 
electronic  charges  of  the  dry  gas  molecules  under  the  influence 
of  an  applied  electromagnetic  field,  the  second  term  the  distor- 
tion for  water  vapor,  and  the  third  term  the  effect  of  the 
orientation  of  the  dielectric  dipoles  of  water  vapor  under  the 
influence  of  a field."  Consistant  with  equation  (3.1)  ki  = 77.6. 
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k2  and  k.3  were  determined  by  Birnbaum  and  Chatterjee  (reference  2), 
who  obtained  k2  = 72.0,  ka  = 3.75  x 10^.  Since  the  total  pressure 
P = Pd  + equation  (3.2)  may  be  rewritten  in  the  form 

N = ki  £ + ka’  I2. 

where  ka'  = ka  - (ki-k2)T.  Since  kj  and  k2  are  nearly  equal,  the 
temperature  variation  of  ka’  is  so  small  that  for  the  region  of 
applicability  of  the  formula  this  may  be  neglected.  Formula  (3.1) 
thus  results. 


3.2  Pressure  and  Temperature  Variation  with  Height 

In  this  section  the  derivation  of  the  required  formulae 
will  closely  follow  that  given  by  Hopfield  (references  3,  4,  5). 
Theoretically,  Ot  , the  rate  at  which  the  temperature  of  dry  air 
decreases  with  height  (the  dry  adiabatic  lapse  rate)  is  given  by 


a,  = ^ SS 


(3.3) 


where 


g = 981  cm/ sec*,  is  the  acceleration  due  to  gravity 


m = 29.0  grams,  is  the  mass  of  a mole  of  dry  air,  and 

R = 8.31  X 10*ergs/ (mole) (deg  Celsius)  is  the  Universal 
Gas  Constant. 

The  theoretical  rate,  9.8°C/km,  is  rather  higher  than  the  observed 
rate,  which  is  closer  to  6°C/km.  The  observed  rate  does,  however, 
fluctuate.  Hopfield  realized  that  if  she  could  choose  a lapse 
rate  a,  which  is  an  integral  fraction  of  (gm/R) , then  her  tropo- 
rpheric  model  would  be  much  simplified.  The  Hopfield  value, 

6.84°  C/km,  corresponds  to 

a «=  I fS.  (3.4) 


The  temperature  lapse  rate  is  assumed  to  be  constant  so  that 


I > 


T = T - ah 
0 


(3.5) 


where  is  the  surface  temperature  and  h the  height  above  the 
surface.  Since,  clearly  the  temperature  cannot  drop  below  absolute 
zero,  the  Hopfield  troposphere  must  possess  an  upper  limit,  hi, 
given  by 


h.  = IJo. 


(3.6) 


Assuming  the  earth  to  be  flat  and  the  acceleration  due  to  gravity 
to  be  constant,  we  deduce  that 


i = - PS. 


(3.7) 


where  p is  the  density  of  air.  Writing  the  gas  equation  PV  = RT 
in  the  form 


= E® 


RT' 


(3.8) 


we  obtain  with  the  aid  of  equation  (3.5) 


i - - gm/R(T„  - ah) . 


whence , 


P = Pq(1  - ^ “/Tq) 


gm 

Ra 


(3.9) 


From  the  above  and  equations  (3.4),  (3.5)  and  (3.6)  we  obtain 


p = p^(l  - h/hi) 


(3.10) 


and 


T = T^(l  - h/hi) 


(3.11) 
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J 


Since  0°C  corresponds  to  273.16°K,  it  follows  that  equation 
(3.6)  may  be  written  as 


hi  = (273.16  + T^,)/a. 

where  is  the  surface  temperature  in  degrees  Celsius.  Using 
the  value  of  a from  equation  (3.4)  we  would  obtain 

hi  = (39.3  + 0.146  T^)km. 

The  actual  formula  employed  by  the  Hop field  model  is 

hi  = (40.1  + 0.149  T^)km.  (3. 


3.3 


The  Zenith  Integral  for  Dry  Air 


The  zenith  integral  is  defined  by 

.hi 


= f N dh, 


(3. 


where  hi  is  the  upper  limit  of  the  modeled  atmosphere.  Since 
e = 0 for  dry  air  it  follows  from  equations  (3.1),  (3.10)  and 
(3.11)  that 

.hi 

= 11^0  I 

'z 


0 


(l-h/hl)“dh, 


i.e.  , 


I = 
z 


77.6  p^hi 

5 T„ 


(3. 


It  is  interesting  to  note  that  if  hi  were  given  by  equations 
(3.6)  and  (3.4),  equation  (3.14)  reduces  to 

= 77.6  p^R/gm,  (3 

so  that  the  zenith  integral  is  a function  solely  of  the  surface 


12) 


13) 


14) 


.14’) 
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pressure.  Assuming  that  the  Earth  is  spherical  and  that  the  inverse 
square  law  of  gravitation  applies,  this  result  may  also  be  deduced 
as  follows.  If  r denotes  the  radial  distance  from  the  Earth's 


center,  then  for  static  equilibrium  the  net  surface  force  on  a 
small  element  of  volume  must  balance  the  body  force.  Hence, 


where  g^  is  the  acceleration  due  to  gravity  at  distance  r^  from 
the  Earth's  center.  Integrating  the  above  equation  between  r^  and 
r,  the  upper  limit  of  the  Earth's  atmosphere,  we  deduce  that 


/: 


Pdr  = p^/g^ 


It  follows  from  equations  (3.1)  and  (3.8)  that  the  'dry'  refractivity 
is  given  by  N = 77.6  pR/m,  whence  the  zenith  integral 


= 77.6  P^R/gom. 


(3.14") 


It  can  thus  be  seen  that  the  result  is  fairly  general  and  is  not 
a function  of  the  temperature  distribution  within  the  atmosphere. 

However,  hi  in  the  Hopfield  model  is  given  by  equation 
(3.12)  rather  than  by  equation  (3.6).  The  corresponding  zenith 
integral  I,  is  given  by 


1 + T^/269 

“ ^'^®Po  1 + T 7273' 

c 

thus,  showing  that  in  the  Hopfield  model  the  zenith  integral 
increases  with  temperature. 


(3.15) 
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0 is  the  location  of  the  receiving  station 
C is  the  center  of  curvature  at  the  point  0, 
X is  the  location  of  the  transmitter, 

OC  •=  To  , the  radius  of  curvature 
QC  *=  r 

OQ  = s,  the  path  length 

angle  QOH  = £„ , the  elevation  angle  at  the  station 
angle  OAC  = %ir. 


r 


The  Ranee  Correction 


The  computed  range  correction  is  given  by 
-ro  +h  1 


As  = 


(n-1)  ^ dr. 


(3.16) 


where  r„  is  the  Earth's  radius  of  curvature  at  the  receiving 
station  and  ds  is  an  element  of  path  length.  The  error  in  performing 
the  integration  along  the  geometric  straight  line  rather  than 
along  the  actual  path  is  small  (reference  3)  "reaching  10%  of  As 
at  the  horizon  but  becoming  neglibly  small  at  elevation  angles 
above  3° or  4° . " 

Referring  to  Figure  3.1,  it  can  easily  be  seen  that 
OA  = ro  sin  E,  and  CA  = rgCOsEo . Since  by  Pythagoras 
QC^  = AQ^  + CA^  it  follows  that 


r^  = (s  + roSinEo)^  + (roCosEo)' 


Hence 


r = (s  + roSinEo);3|-, 


- (r- cosEt 


(3.17) 


Since  N = (n-l)lO®,  it  follows  from  equations  (3.1), 
(3.10),  (3.11),  (3.16)  and  (3.17)  that  the  computed  'dry'  correction 


Asl,  is  given  by 


(l-h/h,)»  r dr 

Asi  = 77.6  X 10" ® ^ 2 -z „ .2  , 

To  J /r^  -(roCOsEo)^ 


(3.18) 


where 


h * r - r„ 


(3.19) 


Writing , 


hi  = r 1 - To  , (3.20) 

equation  (3.18)  may  be  re-written  as 

Asi  = 77.6  X 10-®  ^ f (3.21) 

° r \ 7 /r^-(roCOsE„)2 

O 

As  pointed  out  by  Yionoulis  (reference  6),  equation  (3.21)  may  be 
integrated  in  closed  form,  but  the  numerical  solution,  even  in 
double  precision,  is  inaccurate  at  high  elevation  angles.  Yionoulis* 
answer  to  the  problem  was  to  produce  two  separate  solutions,  one 
being  valid  for  all  elevation  angles  except  the  low  ones,  the  other 
for  all  angles  except  the  high  ones.  A simpler  solution  is 
obtained  by  using  the  NWL  approximation  (Reference  7)  . This  will 
now  be  described. 


Let 


Ni  = 77.6  p„/To  10' 


(3.22) 


Since, 


ri-  r Jri^-  rM/ri+  ro 


ri-  ro  I j. 


(3.23) 


and  since  (ri  + ro)/ri  + r)  is  approximately  equal  to  one,  it 
follows  that  As^  is  approximately  equal  to 

dr 


4S,  - N,  f' 

rn  \^1  - /r^-rc 


(3.24) 


X-Q  \*-  i ~ *-0/  “iO  COS^Eo 

A more  rigorous  justification  for  the  NWL  approximation  (3.24)  is 
given  in  Appendix  B. 

Making  the  substitutions 


X 

= r^  - 

ro^  cos^Eo 

Xi 

= 

- ro^cos^E 

Xo 

- ro^  cos^E 

(3.25) 
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we  deduce  that 


■ 


A closed  form  solution  for  the  above  integral  (reference  7)  is 
easily  obtained.  However,  if  As  i is  computed  by  the  usual  method 
of  evaluating  the  solution  at  both  limits,  then  the  result  will 
be  numerically  inaccurate,  particularly  at  high  elevation  angles. 

A better  way  to  evaluate  the  integral  is  given  in  Appendix  A. 

It  follows  from  equations  (3.26)  and  (3.25),  and  equations  (1),  (3), 
(12)  and  (13)  in  Appendix  A that, 

Asx  = ^ Nilod  + §t  + jt"  + ^t^  + x^f*).  (3.27) 

where 

lo  = (ri-ro) (ri+ro)/d 

d = /(ri-ro)  (ri+ro^)  + ro^  sin^Eo  + rosinEo  ■ (3.28) 

t = lo/d 

since  d^  ^ (ri-ro)  (ri+ro)  , and  ri-ro^O, 

I if  follows  that 

t 

i 

' 0<t;$l  (3.29) 


Since  the  last  term  on  the  right  hand  side  of  equation  (3.27)  con- 
tributes less  than  0.5%  to  the  total  integral,  it  may  be  neglected. 
We  may  therefore  write  As^  in  the  truncated  form 


i.e.  , 

hi  - ho  =[40.1  + 0.6  ho/hi  + 0.149  T,]  km  (3.33) 

Since  hi  is  of  the  order  of  40km,  the  ratio  ho /hi  is 
small  and  may  be  ignored.  Following  O'Toole  (reference  7)  we 
thus  obtain 

hi  - ho  = [40.1  + 0.149  Tc]km,  (3.34) 

In  other  words,  the  height  of  the  atmosphere  above  the 
station  is  independent  of  the  station  height  above  sea-level. 

Taking  station  height  into  account,  equation  (3.21)  must 
be  modified  to 


whence,  by  equations  (3.31)  and  (3.32) 


As.  . 77.6  . 10-‘ 

*^r  o+ho 

Let  So  = To  + ho,  and  ki  = hi  - ho 
Since,  by  equation  (3.20),  hi  = ri 

hi  - ho  = r 


hi**  r dr  /o  35 \ 

(hi-ho) Vr^-(ro+ho)^cos^Eo 

(3.36) 

To,  it  follows  that 
- So.  (3.37) 


We  hence  deduce  that 

AS.  = 77.6  X lO-'IJ  y 

S 0 

Comparing  equations  (3.34),  (3.37)  and  (3.38)  with 
equations  (3.12),  (3.20)  and  (3.21)  we  see  that  the  only  difference 
between  the  equations  is  that  the  radius  of  curvature  tq  has 
been  replaced  by  the  new  radius  of  curvature  Sq.  Since,  further- 
more, the  range  correction  is  comparatively  insensitive  to 
variations  in  the  radius  of  curvature,  we  conclude  that  as  long 
as  the  pressure  and  temperature  are  obtained  at  the  station, 
the  station  height  above  sea- level  may  be  ignored. 


r or 


/r^- (sqCosEo) 


(3.38) 


3 . 6 The  'Wet*  Correction 

The  partial  pressure  of  water  vapor  appearing  on  the 
right  hand  side  of  equation  (3.1)  is  computed  as  a function  of 
the  fractional  relative  humidity  h^.  (h^  = 1 corresponding  to 
100  percent  relative  humidity)  and  the  saturation  water  vapor 
pressure  e . It  is  given  by  the  formula 

e = hj.  6g.  (3.39) 
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where 


eg  = exp(l.  80910  + . (3.40) 

Equation  (3.40)  comes  from  the  National  Climatic  Center  in  Asheville, 
N.C.,  (Hopfield,  Private  Communication,  1976). 

It  was  observed  by  Hopfield  (reference  4)  that  "at 
heights  of  9 or  10km,  the  atmospheric  pressure  is  still  almost 
one  third  of  its  surface  value,  but  the  partial  pressure  of  water 
vapor  is  nearly  zero". 

Seemingly  without  any  theoretical  justification,  Hopfield 
then  went  ahead  to  give  the  'wet'  correction  in  a form  similar 
to  the  'dry'  correction,  equation  (3.12)  being  replaced  by 

ha  = 12.0  km,  (3.41) 

and  equation  (3.22)  by 

Na  = .373  e/To^  , (3.42) 

the  remaining  equations  being  left  the  same  except  for  the 
subscript  1 being  replaced  by  2. 

In  the  original  Hopfield  model,  the  'wet'  atmosphere 
was  assumed  to  have  a ceiling  of  12km  above  sea  level.  In  the 
O'Toole  version  the  ceiling  is  at  12km  above  the  station!  The 
problem  is  really  that  it  is  extremely  difficult  to  modify  a 
model  whose  coefficients  do  not  have  a clear  physical  meaning. 

In  photonap  the  ceiling  is  also  taken  at  12km  above  the  station. 
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3.7 


The  Radius  of  Curvature 


It  is  well  known  (see,  e.g.,  reference  8)  that  the 
radius  of  curvature,  Ra,  at  a point  on  the  Earth's  surface,  must 
satisfy 

ST  " 57  sin^A,  (3.43) 

where 

Rn  = a(l-e^)  (l-e^sin^L)“  » 

Re  = a(l-e^sin^L) 

a = 6378km  is  the  semi-major  axis 

e = 0.081813  is  the  eccentricity 

L is  the  latitude,  and 

A is  the  azimuth. 

It  can  thus  be  seen  that  the  radius  of  curvature  is  a 
maximum  (denoted  by  Rmax ) the  pole  and  a minimum  (denoted  by 
Rj,,,n  ) in  the  North-South  direction  at  the  Equator.  We  find  that 

Rmax  = a(l-e^)“%  = 6400km  (3.44) 

Rmin  = a(l-e^)  = 6335km  (3.45) 

The  average  radius  R is  given  by 

R = a(l-e^)‘  = 6371km  (3.46) 

In  considering  the  Earth's  curvature  we  observe  that  (i) 
it  has  no  effect  on  ray  paths  through  zenith,  and  (ii)  the 
proportional  effect  Increases  as  the  elevation  decreases.  The 
maximum  effect  can  therefore  be  seen  to  occur  at  zero  elevation 
angles.  It  can  be  seen  from  equation  (3.28)  that  for  Eo  = 0, 
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lo  = /ri^-ro^  and  t = 1, 

whence  we  deduce  from  equation  (3.27)  and  equation  (19)  in  Appendix 
A that  the  correction 

Asi  = Ni  /r  -ro^  (3.47) 

Since  ri-ro  = hi  it  follows  that 

Asi(ro)  = N 1 /h  1 (2r o+h 1 ) , (3.48) 

where  we  have  written  the  left  hand  side  in  a form  that  indicates 
that  it  is  a function  of  ro . Differentiating  equation  (3.48) 
logarithmically  we  obtain 


Asi|r7y  3?7  “ 2^0+hi 

From  the  above  we  deduce  that  approximate  relationship 


(3.49) 


6 [As  1 (ro)]  ^ 6rp 
Asi(ro)  2ro+hi 


(3.50) 


choosing  rp  to  be  the  average  radius  R given  by  equation  (3.46), 
we  see  from  equations  (3.44)  and  (3.45)  that  6rp  must  be  numerically 
less  than  36km  for  any  point  on  the  Earth's  surface.  We  hence 
deduce  that  if  we  assume  that  the  radius  of  curvature  is  a constant 
then  the  proportional  error  (the  left  hand  side  of  equation  (3.50)) 
is  less  than  .003.  In  other  words,  taking 


rp  = 6371km 


(3.51) 


leads  to  an  error  that  does  not  exceed  0.3  percent. 
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3.8  Stumnary 

The  following  equations  are  required  for  the  computation 
of  the  range  correction. 


Ni  = 77.6  X 10'®  po/To,  (3.22) 

where  po  is  the  total  atmospheric  pressure 
in  millibars,  and  To  is  the  temperature 
(degrees  Kelvin) 


Na  = .373  e/To^ , (3.42) 

e = h,  exp  f 1.80910  + (3.39) 

where  h,  is  the  fractional  relative  humidity, 
and  Tc  (=To-273. 16)  is  the  temperature 
(degrees  Celsius) 

hi  = 40.1  + 0.149  Tc  (3.12) 

h^  = 12.0  (3.41) 


As  = Asi  + ASa,  is  the  total  correction,  where 

. = N,diti  (^  + l^t,  + l^t,"'  + ^ti^),  for  i = 1,2,  (3.30) 


As 


where. 


di  = /hi(2ro+h,)  + (rosinEo)^  + (rosinEo) 
ti  = hi  (2ro+hi  )/di* 


(3.28) 

(3.20) 


ro  *=  6371,  (3.51) 

and  Eo  is  the  elevation  angle  of  the  transmitter  as  seen  from 
the  station. 

Range  corrections,  computed  for  different  elevation  angles 
and  temperatures,  are  given  in  Table  3.1. 
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TABLE  3.1 


TROPOSPHERIC  RANGE  CORRECTIONS  (HOPFIELD  MODEL) 


Elevation 

Angle 

'Drv' 

correction  (meters) 

'Wet' 

correction  (meters) 

^°C 

-30^C 

^FC 

30^C 

40^C 

-30^C 

^FC 

^0«C 

40“C 

0° 

94.2 

88.3 

83.4 

79.2 

78.0 

0.5 

4.8 

27.3 

44.4 

1° 

63.9 

61.3 

59.0 

57.0 

56.4 

0.3 

2.7 

15.1 

24.6 

2° 

46.6 

45.4 

44.2 

43.2 

42.9 

0.2 

1.8 

9.9 

16.1 

3° 

36.0 

35.3 

34.7 

34.2 

34.0 

0.1 

1.3 

7.2 

11.7 

4° 

29.0 

28.7 

28.3 

28.0 

27.9 

0.1 

1.0 

5.6 

9.1 

6° 

20.6 

20.5 

20.4 

20.3 

20.2 

0.1 

0.7 

3.8 

6.3 

8° 

15.9 

15.9 

15.8 

15.8 

15.8 

0.1 

0.5 

2.9 

4.8 

0 

0 

12.9 

12.9 

12.9 

12.9 

12.9 

0.0 

0.4 

2.4 

3.8 

15° 

8.8 

8.8 

8.8 

8.8 

8.8 

0.0 

0.3 

1.6 

2.6 

0 

0 

CM 

6.7 

6.7 

6.7 

6.7 

6.7 

0.0 

0.2 

1.2 

2.0 

30° 

4.6 

4.6 

4.6 

4.6 

4.6 

0.0 

0.1 

0.8 

1.3 

0 

0 

3.6 

3.6 

3.6 

3.6 

3.6 

0.0 

0.1 

0.6 

1.0 

0 

0 

2.7 

2.7 

2.7 

2.7 

2.7 

0.0 

0.1 

0.5 

0.8 

90° 

2.3 

2.3 

2.3 

2.3 

2.3 

0.0 

0.1 

0.4 

0.7 

The  'dry'  correction  (Asj)  has  been  computed  for  a pressure  of  1013mb. 
For  other  pressures  the  corrections  should  be  adjusted  proportionately. 
The  'wet'  correction  (AS2)  applies  for  100  percent  relative  humidity. 
The  total  correction  (As)  is  given  by 

As  = Asi  + h^As2 

where  h^  is  the  fractional  relative  humidity. 


A. 0 Changes  to  the  Photonap  User’s  Guide 

The  most  important  changes  to  the  User's  Guide  stem  from 
the  added  capability  to  process  geoceiver  measurements.  Geoceiver 
measurements  have  been  assigned  measurement  type  number  27  (see 
Appendix  IV).  Photonap  has  the  capability  to  process  both  'satellite- 
to-ground'  and  ' satellite- to-satellite ' geoceiver  data,  the  User  having 
to  specify  the  type,  which  is  being  processed  (see  category  701  card, 
key  8).  Photonap  computes  ' satellite- to-ground'  tropospheric  corrections 
using  the  NVJL-Hopf ield  model.  The  User  may  be  content  to  use  default 
values  for  meteorological  data  or  he  may  input  his  own  (see  category 
610  card.  The  category  610  card  is  also  used  to  specify  the  geoceiver 
time  interval) . The  geoceiver  time  interval  may  also  be  specified 
as  part  of  the  observation  record  (see  Appendix  I) . Geoceiver  data 
(see  Appendix  IV)  may  be  used  to  recover  orbital  parameters,  station 
locations,  refraction  parameters,  measurement  bias,  measurement  timing 
bias,  measurement  drift  rate  and  measurement  scale. 

Some  relatively  minor  changes  from  NAP  3. IF,  which  give 
Photonap  a Monte  Carlo  capability,  have  been  added  to  Photonap 
(for  more  details,  see  Appendix  B of  the  'Use  of  the  GPS  Satellite 
System  for  the  Determination  of  the  MAGSAT  Position'  NAS5-23587, 

Georg  Morduch,  Old  Dominion  Systems,  September  1976).  Affected 
cards  are  601  cards  (Note  9)  and  604/605  cards. 

The  User  may,  in  the  current  version  of  Photonap,  request 
that  the  covariance  of  the  solution  not  be  computed  (see  Category 
105,  set  2) . 

Description  of  the  mode  2 discrete  thrust  has  been 
clarified  (see  Category  208  card) . 


27- 


am 


References 


1.  Smith,  Ernest  K.  and  Weintraub,  Stanley,  "The  Constants  in 
the  Equation  for  Atmospheric  Refractive  Index  at  Radio 
Frequencies".  Proc.  IRE,  pp.  1035-1037,  August  1953. 

2.  Bimbaum,  G.  and  Chatterjee,  S.K.,  "The  Dielectric  Constant 
of  Water  Vapor  in  the  Microwave  Region",  J.  Appl.  Phys . , 23, 
p.  220,  February  1952. 

3.  Hopfield,  H.S.,  "The  Effect  of  Tropospheric  Refraction  on 
the  Doppler  Shift  of  a Satellite  Signal",  J.  Geophys . Res. 

68,  No.  18,  pp.  5157-5168,  September  15,  1963. 

4.  Hopfield,  H.S.,  "Two-quartic  Tropospheric  Refractivity  Profile 
for  Correcting  Satellite  Data",  J.  Geophys.  Res.,  74,  No.  18, 
pp.  4487-4499,  August  20,  1969. 

5.  Hopfield,  H.S.,  "Tropospheric  Effect  on  Electromagnetically 

Measured  Range:  Prediction  from  Surface  Weather  Data", 

Radio  Science,  6,  No.  3,  pp.  357-367,  March  1971. 

6.  Yionoulis,  S.M.,  "Algorithm  to  Compute  Tropospheric  Refraction 
Effects  on  Range  Measurements",  J.  Geophys.  Res.,  2.5,  No.  36, 
p.  7636,  December  20,  1970. 

7.  O'Toole,  James  W.  , 'Celest  Computer  Program  for  Computing 
Satellite  Orbits'.  Document  NSWC/DL  TR-3565,  Naval  Surface 
Weapons  Center,  Dahlgren,  VA. , October  1976. 

8.  Hosmer,  George  L. , "Geodesy",  John  Wiley  & Son,  Inc.,  New 
York,  1929. 

9.  Morduch,  Georg  E. , "Introduction  to  Orbit  Determination", 

Old  Dominion  Systems,  Inc.,  Final  Report,  Contract  NAS5-20532 
(Task  y/1),  October,  1975. 


-28- 


i 


APPENDIX  A 


EVALUATION  OF  THE  INTEGRAL  I 


n 


The  integral  I^  is  defined  by 


Let 


and 


t = 


■^1-  -^o 


^l+/x^ 


Differentiating  equation  (2)  we  obtain 


1 dx 


From  equation  (3)  we  deduce  that 

2 >6c. 


(1  + t)/t  = 


It  hence  follows  from  equations  (2)  and  (3)  that 
[y  - l][y  - 1 + (l+t)/t]  = 


X - 


(/x^-i/x^)2 


whence, 


^ - [y  - l][l  + yt] 


X - x_. 
X 


Changing  the  variable  of  integration  from  x to  y we  obtain 
the  aid  of  equations  (2),  (4)  and  (7), 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

with 
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1 

Ijj  = y*  (l-y)”  (l+yt)”dy 


Since 


n 


(i^yo"=  E freWt 


(8) 


(9) 


and 


1 

/ 


s=0 


ifS  _ nisi 
^ ^ y (n+s+1) 1 


(10) 


it  follows  that 

I_  = 


n 


TH^T^TTfei+ITT 

s=0 


(11) 


[Note  that  the  definite  integral  of  equation  (10)  is  the  well 
known  Beta  function  and  is  denoted  by  B(n+l,s+l)] 

Evaluating  equation  (11)  for  n = 0,  4 and  5,  we  find  that 


I_  = *6c.  - 
o 1 o 

1 + It  + l-t^  + + 


~ 5^o(^  + + 7^' 

" E^o(^  + yt  + + J^t"  + 5^t*j 

For  the  special  case  t = 1,  equation  (8)  becomes 


(12) 

(13) 

(14) 


j: 


(i-y)”(i+y)”  dy 


(15) 


Since  (1-y) (1+y)  is  an  even  valued  function,  we  may  also  write 
equation  (15)  in  the  form 

/I 

(l-y)"(l+y)"  dy  (16) 
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Making  the  substitution 


u = j(l+y) 


we  obtain 


I = 2^^^  ✓5^,  f u"(l-u)"  du 


The  integral  on  the  right  hand  side  of  equation  (17) 
is  the  form  of  equation  (10).  We  therefore  conclude  that  if 


t = 1. 


In  = 22n(n!)V(2n+l) 


In  particular 


I = /x 


1 

J 
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APPENDIX  B 


THE  NWL  APPROXIMATION  TO  THE  HOPFIELD  INTEGPJ^L 
In  this  appendix  it  will  be  shown  that  equation  (3.24) 


is  a valid  approximation  of  equation  (3.21).  Denoting  the 
difference  by  D,  we  find  that 


whence  G 0 ihd 


D > 0 (5) 


To  obtain  an  upper  bound  of  D,  we  write  G in  the  form 


I 


Hence , 


G < 4 


(Sft)  [‘  ■ Si-] 


(7) 


It  hence  follows  for  equations  (2)  and  (7)  that 

•'ro 


r dr 

/ri^-r^\  /ri^-rM 

/r^-ro^  cos^Eo 

\ri^ -ro7 

(8) 


Making  the  substitutions 

X = - ro^cos^Eo,  Xo  = r -ro^  cos^Eq  , Xi  = r -ro^  cos^Eo , 

inequality  (8)  may  be  written  in  the  form 


(9) 


where 


_ 1 f ^ [~xi  -X  1^ 

Xo 


(10) 


Since  the  total  correction  As i may  be  written  in  the  form 


ASi  = Ni  I4, 


(11) 


it  follows  that  the  proportional  error  in  the  correction  6si,  (=D/Asi 
must  satisfy 


6s  1 <4 


From  equations  (13)  and  (14)  in  Appendix  A,  we  obtain 

I5  5(1  + ait  + aat^+  aot^  + a4t“  + ast®  ) 

IT  “ 6a  + bit  + b2t^+  bTt'  + b.it'*  + bst^  T 


(12) 


(13) 
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where  ai  = 5/7,  az  = 5/14,  as  = 5/42,  ai,  = 1/42,  as  - llh62. 


bi  = 2/3,  bz  = 2/7,  bs  = 1/14,  b4  = 1/126,  bs  = 0. 

Since  for  each  s,  a > b , it  follows  that  the  ratio  Is/I-,  is 

s s 

an  increasing  function  of  t.  Consequently,  Is/Iu^  5/6,  and 


i 


e.  , 


6s  I 


< 


i 


ri-rp 

r i-ro 


6s  1 


2(ri-ro) 

3(ri+ro) 


(14) 


Since  rj  > rp  , we  finally  obtain 

6si  < (15) 

Using  the  values  ri-rp  = 50km,  rp  = 6000km,  we  conclude  that  the 
NWL  approximation  to  the  Hopfield  formula  is  accurate  to  at  least 
one  third  of  one  percent. 
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